{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.3 Discontinuous Galerkin Discretizations for linear transport\n",
    "We are solving the scalar linear transport problem\n",
    "\n",
    "Find $u: [0,T] \\to V_D := \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = u_D\\}$, s.t.\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega} \\partial_t u v +  b \\cdot \\nabla u v = \\int_{\\Omega} f v \\qquad \\forall v \\in V_0 = \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = 0\\}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen import gui"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As a first example, we consider the unit square $(0,1)^2$ and the advection velocity $b = (1,2)$. \n",
    "Accordingly the inflow boundary is $\\Gamma_{in} = \\{ x \\cdot y = 0\\}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.geom2d import SplineGeometry\n",
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (1, 1), \n",
    "                 bcs = (\"bottom\", \"right\", \"top\", \"left\"))\n",
    "mesh = Mesh( geo.GenerateMesh(maxh=0.125))\n",
    "Draw(mesh)\n",
    "order = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convection and boundary values:\n",
    "\n",
    "* To make cases with `CoefficientFunctions` we use the `IfPos`-`CoefficientFunction`. `IfPos(a,b,c)`: `a` decides on the evaluation. If `a` is positive `b` is evaluated, otherwise `c`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "b = CoefficientFunction((1+sin(4*pi*y),2))\n",
    "Draw(b,mesh,\"wind\")\n",
    "from ngsolve.internal import visoptions; visoptions.scalfunction = \"wind:0\"\n",
    "ubnd = IfPos(x-0.125,IfPos(0.625-x,1+cos(8*pi*x),0),0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We consider an Upwind DG discretization (in space):\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v - b \\cdot \\nabla v u + \\int_{\\partial T} b_n \\hat{u} v = \\int_{\\Omega} f v, \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\hat{u}$ is the Upwind flux, i.e. $\\hat{u} = u$ on the outflow boundary part $\\partial T_{out} = \\{ x\\in \\partial T \\mid b(x) \\cdot n_T(x) \\geq 0 \\}$ of $T$ and $\\hat{u} = u^{other}$ else, with $u^{other}$ the value from the neighboring element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "There is quite a difference in the computational costs (compared to a standard DG formulation) depending on the question if the solution of linear systems is involved or only operator evaluations (explicit method). \n",
    "\n",
    " * **We treat the explicit case here only**\n",
    " * and refer to [unit 2.8](../unit-2.8-DG/DG.ipynb) for solving linear systems with DG formulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Explicit time stepping with a DG formulation\n",
    "\n",
    "Explicit Euler:\n",
    "$$\n",
    "\\sum_{T} \\int_T u^{n+1} v = \\sum_{T} \\int_T u^{n} v + \\Delta t \\sum_{T} \\left\\{ \\int_T  b \\cdot \\nabla v u \n",
    "+ \\int_{\\partial T} b_n \\hat{u} v \\right\\} + \\Delta t \\int_{\\Omega} f v, \\quad \\forall v \\in V_h,\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "  M u^{n+1} = M u^{n} - \\Delta t C u^n + \\Delta t f\n",
    "$$\n",
    "\n",
    "In our first example we set $u_0 = f = 0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing convection applications $C u^n$\n",
    "* We can define the bilinear form **without setting up a matrix** with storage. (`nonassemble=True`)\n",
    "* A `BilinearForm` is allowed to be **nonlinear in the 1st argument** (for operator applications)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VT = L2(mesh,order=order)\n",
    "u,v = VT.TnT()\n",
    "c = BilinearForm(VT, nonassemble=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* To access the **neighbor** function we can use `u.Other()`, cf. [unit 2.8](../unit-2.8-DG/DG.ipynb).\n",
    "* To incorporate **boundary conditions** ($\\hat{u}$ on inflow boundaries), we can use the argument `bnd` of `.Other()` (not possible in [unit 2.8](../unit-2.8-DG/DG.ipynb) )\n",
    "* If there is no neighbor element (boundary!) the `CoefficientFunction` `bnd` is evaluated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "n = specialcf.normal(mesh.dim)\n",
    "upw_flux = b*n * IfPos(b*n, u, u.Other(bnd=ubnd))\n",
    "dS = dx(element_boundary=True)\n",
    "c += - b * grad(v) * u * dx + upw_flux * v * dS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bilinearform for operator applications\n",
    "Due to the `nonassemble=True` in the constructor of the `BilinearForm` we can use `c.mat * x` to realize the operator application without assembling a matrix first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(VT); \n",
    "Draw(gfu,mesh,\"u\",sd=4,min=-0.1,max=2.1,autoscale=False)\n",
    "res = gfu.vec.CreateVector()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Operator application (equivalent to assemble and mult but faster)\n",
    "res.data = c.mat * gfu.vec "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# this does not work:\n",
    "c2 = BilinearForm(VT)\n",
    "res.data = c2.mat * gfu.vec "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving mass matrix problems (see also [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb))\n",
    "\n",
    "* Need to invert the mass matrix\n",
    "* For DG methods the mass matrix is block diagonal, often even diagonal.\n",
    "* FESpace offers (if available):\n",
    "  * `Mass(rho)`: mass matrix as an operator  (not a sparse matrix)\n",
    "  * `Mass(rho).Inverse()`: inverse mass matrix as an operator (not a sparse matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "invm = VT.Mass(1).Inverse()\n",
    "print(invm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The (simple) time loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "#better: %%timeit\n",
    "t = 0; gfu.vec[:] = 0\n",
    "dt = 2e-4\n",
    "tend = 0.6\n",
    "\n",
    "while t < tend-0.5*dt:\n",
    "    res.data = invm @ c.mat * gfu.vec     \n",
    "    gfu.vec.data -= dt * res\n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Skeleton formulation (cf. [unit 2.8](../unit-2.8-DG/DG.ipynb))\n",
    "\n",
    "So far we considered the DG formulation with integrals on the boundary of each element. Instead one could formulate the problem in terms of facet integrals where every facet appears only once. \n",
    "\n",
    "Then, the corresponding formulation is:\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v - b \\cdot \\nabla v u + \\sum_{F\\in\\mathcal{F}^{inner}} \\int_{F} b_n \\hat{u} (v - v^{neighbor}) + \\\\\n",
    "  \\sum_{F\\in\\mathcal{F}^{bound}} \\int_{F} b_n \\hat{u} v = \\int_{\\Omega} f v , \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\mathcal{F}^{inner}$ is the set of interior facets and $\\mathcal{F}^{bound}$ is the set of boundary facets. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### integral types:\n",
    "The facet integrals are divided into inner facets and boundary facets.\n",
    "To obtain these integrals we combine the `DifferentialSymbol` `dx` or `ds` with `skeleton=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dskel_inner  = dx(skeleton=True)\n",
    "dskel_bound  = ds(skeleton=True)\n",
    "# or:\n",
    "dskel_inflow = ds(skeleton=True, definedon=mesh.Boundaries(\"left|bottom\"))\n",
    "dskel_outflow = ds(skeleton=True, definedon=mesh.Boundaries(\"right|top\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Skeleton formulation of $C$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = BilinearForm(VT,nonassemble=True)\n",
    "c += -b * grad(v) * u * dx\n",
    "c += upw_flux * (v-v.Other()) * dskel_inner\n",
    "\n",
    "c += b*n * IfPos(b*n,u,ubnd) * v * dskel_bound\n",
    "#alternatively (if you know in/out beforehand):\n",
    "#c += b*n * ubnd * v * dskel_inflow\n",
    "#c += b*n * u * v * dskel_outflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The (simple) time loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "t = 0; gfu.vec[:]=0\n",
    "while t < tend-0.5*dt:\n",
    "    gfu.vec.data -= dt * invm @ c.mat * gfu.vec \n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## TraceOperator (Purpose: avoid neighbor access `Other()`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The trace operator maps from an `FESpace` to a compatible `FacetFESpace` by evaluating the traces. \n",
    "\n",
    "The **stored value** on a facet is the **sum of the traces** of the aligned elements. \n",
    "\n",
    "No need to access the neighbor (`.Other()`) later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VF = FacetFESpace(mesh, order=order)\n",
    "V = FESpace( [VT, VF] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = VT.TraceOperator(VF, False)\n",
    "print(\"The space VT has {} dofs.\".format(VT.ndof))\n",
    "print(\"The trace space VF has {} dofs.\".format(VF.ndof))\n",
    "print(\"trace represents an {} x {} operator\".format(trace.height, trace.width))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gfuF = GridFunction(VF)\n",
    "gfuF.vec.data = trace * gfu.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To make use of the trace operator, we split the previous (`element_boundary`)formulation into three parts:\n",
    "$$\n",
    "  \\sum_{T} \\int_T - b \\cdot \\nabla v u + \\sum_{T} \\int_{\\partial T} b_n \\hat{u}^* v\n",
    "  + \\sum_{F\\in\\mathcal{F}^{inflow}} \\int_{F} b_n u^{inflow} v \\quad \\text{ with } \\quad \\hat{u}^* = f(u^{left}+u^{right},u)\n",
    "$$\n",
    "\n",
    "* part 1: element local volume integrals (no interaction with neighbor elements) (on `VT`)\n",
    "* part 2: couplings on the facets (on `VT` $\\times$ `V`)\n",
    "* part 3: (rhs + ) boundary terms (on `VT`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "uT, vT = VT.TnT()\n",
    "c1 = BilinearForm(space=VT, nonassemble=True)                     # part 1\n",
    "c1 += -uT * b * grad(vT) * dx\n",
    "uT,uF = V.TrialFunction() \n",
    "c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True)   # part 2\n",
    "# here uf-u = u_me + u_other - u_me is the neighbor value (0 at bnd)\n",
    "c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)\n",
    "rhs = LinearForm(VT)                                              # part 3\n",
    "rhs += b*n * IfPos(b*n, 0, ubnd) * vT * dskel_bound # dskel_inflow\n",
    "rhs.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that `c1` and `c2` act on different spaces. To make them compatible, we combine:\n",
    " * `TraceOperator`: VT $\\to$ VF\n",
    " * `Embedding`s: VT $\\to$ VT $\\times$ VF and VF $\\to$ VT $\\times$ VF (cf. [unit 2.10](../unit-2.10-dualbasis/dualbasis.ipynb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The Embeddings are (with $n_F = \\operatorname{dim}(V_F)$, $n_T = \\operatorname{dim}(V_T)$, $n_V=n_F+n_T$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "$$\n",
    "\\texttt{embV}: V_T \\to V_T \\times V_F, \\qquad \n",
    "\\texttt{embV} = \n",
    "\\underbrace{\\left( \\begin{array}{c} I \\\\ 0 \\end{array} \\right)}_{\\texttt{embT} \\in \\mathbb{R}^{n_V \\times n_T}}\n",
    "+ \n",
    "\\underbrace{\\left( \\begin{array}{c} 0 \\\\ I \\end{array} \\right)}_{\\texttt{embF} \\in \\mathbb{R}^{n_V \\times n_F}} \\underbrace{\\operatorname{tr}}_{\\texttt{trace}\\in \\mathbb{R}^{n_F \\times n_T}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "embT = Embedding(V.ndof, V.Range(0))\n",
    "embF = Embedding(V.ndof, V.Range(1))\n",
    "embV = embT + embF @ trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The operators can now be combined easily. Recall:\n",
    " * `c1`: `VT` $\\to$ `VT`' $\\quad$ $\\bullet$ `c2`: `VT` $\\times$ `VF` $\\to$ `VT`' $\\quad$ $\\bullet$ $\\Rightarrow$ `c`: `VT` $\\to$ `VT`'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = c1.mat + c2.mat @ embV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "t = 0; gfu.vec[:] = 0\n",
    "invm = VT.Mass(1).Inverse()\n",
    "r = gfu.vec.CreateVector(); r.data = invm * rhs.vec # <- the boundary data\n",
    "while t < tend-0.5*dt:\n",
    "    gfu.vec.data -= dt * (invm @ c * gfu.vec + r)\n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Careful with one-liners!**\n",
    "\n",
    "The following will not work correctly: \n",
    "(`gfu` appears left and right: order of ops is important!)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "t = 0; gfu.vec[:] = 0\n",
    "while t < tend-0.5*dt:\n",
    "    gfu.vec.data -= dt * (r + invm @ c * gfu.vec)\n",
    "    t += dt; Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In doubt? Create a temporary vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tmp = gfu.vec.CreateVector()\n",
    "t = 0; gfu.vec[:] = 0\n",
    "while t < tend-0.5*dt:\n",
    "    tmp.data = r + invm @ c * gfu.vec\n",
    "    gfu.vec.data -= dt * tmp\n",
    "    t += dt; Redraw(blocking=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Speed up with `geom_free` integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "With $u(x) = \\hat{u}(\\hat{x})$, $v(x) = \\hat{v}(\\hat{x})$ (directly mapped) and $\\mathbf{w}(x) = \\frac{1}{|\\operatorname{det}(F)|} F \\cdot \\hat{\\mathbf{w}}(\\hat{x})$ (Piola mapped) there holds\n",
    "$$\n",
    " \\int_T \\mathbf{w} \\cdot \\nabla u \\ v  = \\int_{\\hat{T}} \\hat{\\mathbf{w}} \\cdot \\hat{\\nabla} \\hat{u} \\ \\hat{v}\n",
    "$$\n",
    "\n",
    "and similarly for the facet integrals. \n",
    "\n",
    "Integrals are independent of the \"physical element\" (up to ordering) similar to integrals in [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb).\n",
    "\n",
    "$\\leadsto$ allows to prepare intermediate computations for a large class of elements at once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Piola wind:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfwind = GridFunction(HDiv(mesh, order=order))\n",
    "gfwind.Set(b)\n",
    "b=gfwind; # <- overwrite old name \"b\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Same as before, but with  `geom_free=True` flag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# part 1\n",
    "uT, vT = VT.TnT()\n",
    "c1 = BilinearForm(space=VT, nonassemble=True, geom_free=True)\n",
    "c1 += -uT * b * grad(vT) * dx\n",
    "# part 2\n",
    "uT,uF = V.TrialFunction() \n",
    "c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True, geom_free=True)\n",
    "# here uf-u = u_me + u_other - u_me is the neighbor value\n",
    "c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = c1.mat + c2.mat @ embV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "t = 0; gfu.vec[:] = 0\n",
    "while t < tend-0.5*dt:\n",
    "    gfu.vec.data -= dt * (invm @ c * gfu.vec + r)\n",
    "    t += dt\n",
    "    Redraw(blocking=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
